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1. Introduction 

A force between nucleons, the nuclear force, is the most fundamental quantity 
in nuclear physics. The phenomenological nucleon-nucleon (NN) potentials^ ex- 
hibit long-to-medium distance attractions, which have been explained by meson 
exchanges between nucleons^ while they have short distance repulsion, the repul- 
sive coreP essential for the stability of nuclei against collapse, whose origin has 
not yet been well understood. Moreover, general short distance repulsions among 
baryons, not only NN but also three- nucleon (3N), baryon-baryon (BB) and three- 
baryon (3B)J31may be required to explain a recently observecP neutron star as heavy 
as twice a solar mass . Unfortunately it is difficult or almost impossible to deter- 
mine these short distance repulsions experimentally, except for the NN case, while 
the theoretical determination for these short distance repulsions requires difficult 
non-perturbative calculations in QCD. 

Recently, a novel method to define the NN potential in QCD has been pro- 
posed*^ (we call it the HAL QCD scheme or method in this paper), and the 
method has been widely applied to various hadronic interactions^^] including BB 
potentials"! and 3N forces:- 7 Furthermore, in the HAL QCD scheme, the short 
distance behavior of two- and three baryon potentials can be investigated analyt- 
ically via the operator product expansion combined with a renormalization group 
(RG) analysis of perturbative QCDpHU] 

In this paper, we review both analytical investigations and numerical results for 
short distance behaviors of these potentials in the HAL QCD method . In Sec. [2] we 
briefly explain the HAL QCD method to define the potential in QCD, taking the NN 
case as an example. The OPE (operator product expansion) formalism to analyze 
the short distance behaviors of potentials in the HAL QCD scheme is summarized 
in Sec. |3l Results for potentials in lattice QCD simulations are summarized and 
compared with the OPE analyses in Sec. [4) and a summary is given in Sec. 

2. HAL QCD method 

In this section, we explain the HAL QCD method to define the potential in QCD 
for the NN case. The basic quantity is the equal-time Nambu-Bethe-Salpeter(NBS) 
wave function, defined by 

P k a T S2 W = (0\T {N a (x o ,O)N (x Q + x,0)}\NN,k, Sl s 2 ) in , (1) 

where (0| is the QCD vacuum (bra)-state, \NN, k, siS2)i n is the two-nucleon asymp- 
totic in-state with helicity si, s 2 and relative momentum k in the center of mass 

frame, whose total energy is given by Wk = 2 \J k 2 + m 2 N with the nucleon mass 
mjv, and T represents the time-ordered product. The local interpolating operator 
for the nucleon is taken as N a (x) = e a bc{u a (x) T C'j^d 1 ' {x))q%{x) (with x = (x,t)), 
where a, b, c are the color indices, a is the spinor index, C = 7274 is the charge 
conjugation matrix, and q(x) = (u(x),d(x)). 
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Short Distance Repulsion Among Baryons 3 

As long as the total energy Wk lies below the pion production threshold(i.e. 
Wk < Wth — 2rriN + with pion mass m w ), the NBS wave function at r = \x\ —> 
oo satisfies the Helmholtz equation: 

[fc 2 + V 2 ]^<»^0, fc = |fe|, (2) 

where c = S1S2 and T = a/3. Moreover, the radial part of the NBS wave function 
for given orbital angular momentum L and total spin S for asymptotically large r 
becomeiPI^U 

V \r;LS) oc Hfa-^ + W^w (3) 
kr 

where <5ls(^) is the NN phase shift obtained from the S-matrix in QCD below the 
inelastic threshold. 

From the above property of the NBS wave function, it is possible to define a 
non-local potential through the equation 

■„2 y2^ 
"27, 



[E k -H ]$ c (x) =J2j U r ., ri (x,y)^l c (y)d 3 y, ^E k = H 



(4) 



with the reduced mass /1 = mjv/2, where U{x,y) is expected to be short-ranged 
because of absence of massless particle exchanges between two nucleons. The po- 
tential U(x, y) is finite and does not depend on the particular renormalization 
scheme, since the NBS wave function cp a 'p is multiplicatively renormalized. While 
Lorentz covariance is lost by using the equal-time NBS wave function and Eq. (j4|) is 
written as a Schrodinger type equation, a non-relativistic " approximation" is NOT 
introduced to define U(x,y). 

The non-local but energy(fc)-independent potential U(x,y) can be formally 
given byP"^ 

U ri , ra (x,y) J £ [^,-^o]^^ Cl N^ Ci;fc2 ^ 2 ^' C2 ( y )t, (5) 

fcl,Cife 2 ,C2 



where fc t h is the threshold momentum defined to satisfy 2 */ fc 2 h + m 2 N = Wth so that 
the summation over fei, is restricted below the inelastic threshold, and N^ 1 is 
the inverse of Af defined from the inner product of NBS wave functions as 

A/ -fc 1 , cl ;fe 2 , C2 = (y=i,c 1)¥ ,fe 2 ,c 3 ) =J2fd 3 x ifft'" 1 ^ <Pr ,C2 (x) (6) 

r •* 

for |fei I , |&2 1 < fc t h- This U(x, y) is energy-independent by construction and satisfies 
eq. (j4|) as 

J2 J U r .,rAx,y)<Pr[ c (y)d 3 y = ''^ " ^ - H ] ^ (x)N£ cyM ^ ■ M k ^ k c 

Ti k 1 ,c 1 k 2 ,C2 

= [E k - H ] 4' c (x) (7) 



May 1, 2013 8:34 WSPC/INSTRUCTION FILE RepulsiveCore3-5 



4 S. Aoki, J. Balog, T. Doi, T. Inoue, P. Weisz 

as long as Wk < Wth- Once the non-local potential U(x,y) which satisfies eq. ((H) is 
obtained, eq. (|3]) guarantees that we can extract the phase shift SLs(k) at Wk < W t h 
in QCD by solving the Schrodinger equation with this potential. Note however that 
the non-local potential which satisfies eq. Q at Wk < W t h is not unique, since 
the potential itself is not a physical observable. We may add terms orthogonal to 
all NBS wave functions below the inelastic threshold, which only affect eq. Q at 
W k > W th ® 

The construction of U(x, y) in eq. ([5]) is formal and is inadequate for practical 
use, since only a limited number of wave functions at low energies (ground state and 
possibly a few low-lying excited states) can be obtained in lattice QCD simulations 
in a finite box. In practice, we therefore expand the non-local potential in terms of 
the velocity (derivative) with local coefficient functions; 2 ^ U(x, y) = V(x, V)S 3 (x — 
y), whose leading order terms reads 

V LO (x) - V (r) + V a (r)a 1 ■ <r 2 + V T {r)S 12 , (8) 

where cri is the Pauli-matrix acting on the spin index of the z-th nucleon, and 
S\2 = 3(x ■ (T\){x ■ <r 2 )/r 2 — cri ■ <r 2 is the tensor operator. It has been found 
in numerical simulations^ that contributions from higher order terms are much 
smaller than those from V LO at low energy. This means that non-locality is rather 
small at low energy in our definition of the potential. 
For example for L — and S = 0, we obtain 

V c (r, L = S = 0) — V (r) - W a (r) = - - - - Q) . (9) 

As mentioned before, since the potential itself is not a physical observable, 
its short-distance behavior depends on its definition, in particular, the choice of 
nucleon operator in the NBS wave function. In this review, we analyze the short 
distance behaviors of the potentials defined in the HAL QCD scheme exclusively, 
by comparing analytic results to numerical results from lattice QCD simulations. 
The method in this review, however, can be easily applied to other definitions for 
the potential. 

3. OPE and short distance behavior of the potential 
3.1. Notation and Method 

One of the outstanding properties of QCD is asymptotic freedom; renormalized 
running couplings g 2 (r) which characterize interactions at short distance scales r 
tend logarithmically to zero as the scale decreases g 2 (r) ~ — 2 ^ o ^n^] , where f3 = 
[11 — 2Nf/3] /(167r 2 ) for gauge group SU(3), provided the number of dynamical 
quarks is limited (Nf < 16). This has the important consequence that the behavior 
of a wide class of physical quantities at short distances (alternatively high energies) 
can be computed perturbatively. The NBS wave functions defined in eq. ([T]) belong 
to this class. In the framework of perturbation theory of QCD it is convenient to 
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regularize the ultra-violet divergences using dimensional regularization (DR), and 
then bare parameters and operators must be renormalized. 

The operator product expansion (OPE) method is based on the observation that 
the product of the two renormalized nucleon operators appearing in the definition 
of the NBS wave function is for short distances expressed as a sum of local gauge 
invariant operators times coefficient functions: 

N a (0, 0)N/,(x, 0) x ~ q Y, F^(x)O a , (10) 

A 

which holds when sandwiched between asymptotic states or in correlation functions 
with fundamental fields at widely separated distances. The coefficient functions 
F^o(x) appearing above can be computed perturbatively using the RG equations 
which express independence of the bare theory on the renormalization scale. □ If the 
nucleon operators and the Oa are chosen so that they renormalize multiplicatively 
without mixing: Oa = ZaOa^htc , then the corresponding coefficients behave as 
^ap( x ) ~ g 2 ^ A ~^ A H\ x \) i where £a determines the loop order at which the corre- 
sponding operator first appears in perturbation theory (PT) and Pa is related to 
the difference of the 1-loop anomalous dimension ja of the operator Oa and that 
of two nucleons 7^ = 1/An 2 : p A = {jA — 27jv)/(2/3o)- 

The leading asymptotic short distance behavior of the NBS wave function is 
then dominated by the operator with the smallest value of Ia — Pa, provided of 
course that its matrix element between the vacuum and the particular 2N state does 
not vanish (which is generically the case). Usually the smallest value of I a — Pa is 
attained for operators Oa which appear at tree level [l A = 0) i.e. c A ^ in the 
free field expansion: 

iV a;bare (0,0)iV^ 0) ~ go= o : NaibareNftbare = (0, 0) + O(a) 

= 5> s O B;baro (0,0) + O( : z). (11) 

B 

Given an arbitrary basis of gauge invariant r-quark operators the anomalous 
dimensions are related to the eigenvalues of the matrix Z specifying the mixing 
under renormalization. The renormalization constant matrix can be determined by 
computing the correlation function of the operators multiplied by r quark fields 
at widely separated distances. The 1-loop computation just involves the diagrams 
with a gluon line connecting a pair of quark lines emanating from the operator 
vertex. The divergent part, using DR in D = 4 — 2e dimensions is symbolically of 
the form: 



, -1 1 — loop,div fj2 1 

(z)$ 9 (z)\ -^7[( T o + ATi)-g a (x)®g b ( a; )]^ (12) 



where To , Ti are matrices in flavor and spinor space. The term proportional to 
the covariant gauge parameter A cancels with the renormalization of the external 

a For details we refer the reader to sect. 3 of refPS 
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quark lines. Thus only To is relevant here, which in a chiral representation of the 
gamma matrices, i.e. 75 = diag(l, 1, —1, —1) , is simply given by 

(13) 

where either € {1, 2} (right-handed) or a\,P\ <E {3, 4} (left-handed). 

Now any local gauge invariant 6-quark operator of (lowest) canonical dimension 
9 can be written as a linear combination of operators of the form 

B™{x) = B*l(x)Bg(x) , B*(x) = flg( S ) = e abc q a a 'f(x) q ^(x)q c y ' h (x) , 

(14) 

where T — {a,/3,j} are sets of spinor labels and F = {f,g,h} are sets of flavor 
labels. 

Due to the structure of (fl3"j) , which reflects the chirality conservation of massless 
QCD, the operators in eq. (|14[) mix only with operators which preserve the set of 
flavor and Dirac indices in the chiral basis i.e. 

f 1 uf 2 = f[ u f' 2 , ri u r 2 = ri u r 2 . 

We thus start by constructing sets of operators with fixed flavor and Dirac struc- 
ture. Note however that such operators are not all linearly independent. Relations 
between them of the form 

3^S+E^? a a=0, (15) 

follow from an identity satisfied by the totally antisymmetric e symbol and the 
Grassmannian nature of the quark fields. Here the z-th index of abc and j-th index of 
def are interchanged in (abc, def)[i,j]. For example, (ri, ra)[l, 1] = o:2/3i7i, cti/?272 
or (ri,r2)[2,l] = aia27i, Pi02l2- Note that the interchange of indices occurs si- 
multaneously for both. 

For a given set of Dirac and flavor indices the initial set of operators may be 
quite large. The determination of an independent basis and the diagonalization 
of the renormalization constant matrix Z is then more conveniently done using 
algebraic computer programs written in Mathematica or Maple. 

Although the NBS wave function is always dominated by the operator with the 
largest value of f3 = fie, for the potential this is only the case if /3c ^ 0. The 
potential is repulsive at short distances if /3c < and attractive if /3c > 0. However 
if /3c — the leading term in the potential is governed by the next largest non- 
zero anomalous dimension, and the sign of the potential requires non-perturbative 
additional information. That there are operators with (3c = is clear; e.g. operators 
of the form B 'J^j a ]B 9 fA - where a, f3 have the same chirality and a, (3 the opposite 
chirality, since there is no contribution from diagrams where the gluon line joins 
quarks in the different baryon operators. 
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We have carried out the program above for the case of 2-baryon wave functions 
involving 3 different flavors of quarks and also for 3-baryon wave functions. The 
main results from these analyses are presented in the next subsection. 

3.2. Results 

We have considered the following four problems: short-distance interaction between 
two nucleonsj^ between two octet baryonspSl among three nucleon^S and finally 
among three octet baryons!^ In all these cases we mainly concentrated on the ques- 
tion of the existence of a repulsive (or attractive) core of the interaction potential. 

3.2.1. 2- body forces 

The technical part of the computation^ for these cases consists of three main steps: 
the enumeration of independent (taking into account the identities fj 1 5|) ) local gauge 
invariant 6-quark operators, calculation of the spectrum of 1-loop anomalous dimen- 
sions (by diagonalizing the 1-loop mixing matrix) and finally checking whether a 
given operator appears already at tree level in the OPE. We found that all 1-loop 
anomalous dimension eigenvalues are, when multiplied by 487r 2 , even integers. This 
suggests that there is a rationally related natural basis in which the mixing matrix 
is automatically diagonal but we were unable to find a simple method to find this 
basis. 

2-nucleon potential (2 flavors}^ 

The most important property of the spectrum in this case is that for all operators 

1A < 2jn (16) 

hence /3a, the leading power of the running coupling is non-positive. Almost all 
operators have negative Pa indicating repulsive behavior, but as explained above, 
there are also operators with the equality sign in (fT6|). In these cases whether the 
leading behavior corresponds to attraction or repulsion cannot unfortunately be 
decided using perturbation theory alone. In the short distance limit the potential 
behaves as 

V(x)~r£ *™ (17) 

ar 

in these cases, where f3' < corresponds to the subleading operator and R is the 
ratio of the matrix elements of the leading and subleading operators sandwiched 
between the vacuum and the given 2-nucleon state. The latter can only be calculated 
non-perturbatively. Since we are only interested in the sign of this quantity, we were 
able to use a simple effective model to conclude^ that this sign is always positive 
leading to repulsion in all cases. 

2-baryon potential (3 flavors}^ 
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Table 1. List of channels with 1-loop anomalous dimensions greater than 2fpf, others can be 
obtained by symmetry transforms of the Dirac indices 1 f> 2,3f> 4 or (1,2) f> (3,4). 



Dirac structure 


48tt 2 7a 


SU(3) representation 


112334 


42 


1 


111222 


36 


1 


112234 


36 


1 


111234 


32 


1®8 


112234 


32 


1®8 



The tensor product of two octets can be decomposed under SU(3) as 

8® 8 = l s ®8 s ffl27 s ffi8 a ffll0 a fflT0 Q , (18) 

where the subscript s(a) indicates symmetry (antisymmetry) under the interchange 
of the two octet representations. We considered local 6-quark operators with flavor 
content uuddss because such operators represent all multiplets in the above decom- 
position. Again, most operators have ja < but there are also attractive cases, 
operators with ja > "^In (see Table QJ. These all belong to the representations l s , 
8 S and 8 a , consistent with what we found for the 2-nuclcon case, since 2-nucleon 
states belong to the 10 a and 27 s SU(3) representations. Finally we found that all 
operators with 7^4 > 277V appear already at tree level in the OPE of two nucleon 
fields. 

3.2.2. 3-body forces 

As explained in the introduction, collective phenomena in nuclear physics cannot 
be described using only 2-body potentials. In particular, the existence of massive 
neutron stars suggests that also the 3-nucleon potential has a repulsive core. To 
study 3-body forces with our method we introduce the 3-particle generalization of 
the NBS wave function ([T]) 

<fgp y (r,p) = (0\T{B a ( Xl ,0)B p {x 2 ,0)B^x 3 ,0)}\W,E) in , (19) 

where for transparency we suppressed momentum and helicity dependence. Here 
B a is the baryon-octet generalization of the nucleon doublet field appearing in 
([I]) and 3B is an appropriate 3-baryon state with energy E. In the argument of 
the 3-body NBS wave function we use Jacobi coordinates r = X\ — x 2 and p = 
[x 3 - {x 1 +x 2 )/2]/V3. 

In the local potential approximation we have 

<3 7 = E^ M . (20) 



m N V 4 "J 



May 1, 2013 8:34 WSPC/INSTRUCTION FILE RepulsiveCore3-5 



Short Distance Repulsion Among Baryons 9 

This defines the total 3-baryon potential V3B and the true 3-baryon potential V3B 
is given by the potential excess 

V 3B (r,p) =^2V aB {xi -Xj) + V 3B (r,p), (21) 

where V2B is the 2-baryon potential discussed above. 

In the UV limit, i.e. when s — \Jr 2 + p 2 — > we can use the OPE method to 
calculate the leading behavior of the wave function (ITO1) analogously to the 2-body 
case but here we have to consider local gauge invariant 9-quark operators. We have 
carried out the anomalous dimension analysis both in the 3-nucleon case and the 
most demanding 3-baryon case. 

3-nucleon potential ( 2 flavors}^ 

The most striking feature of the spectrum in this case is that for all operators 

1A < 3 7 w. (22) 

Here again operators with the largest anomalous dimension are present in the 3- 
nucleon OPE at lowest order and we conclude that the total local potential in this 
case behaves asymptotically as 

VMr.p)-'^- (23) 

(up to a positive constant). Since at short distances this dominates the 2-nucleon 
potential which behaves as given in (| 1 T[) it follows that the genuine 3-body potential 
V3B introduced above is responsible for the short-distance repulsion. 

We note that showing the presence of a repulsive core in the 3-nucleon case is 
the most important result of our perturbative considerations. Universal repulsion 
follows unambiguously from (|22[) without having to rely on any model-dependent 
results. 

We also note that one can observe from the calculated pattern of the spectrum 
of 1-loop anomalous dimensions that the Pauli exclusion principle is at work here: 
increasing the number of quarks while keeping the number of flavors fixed makes 
repulsion stronger while increasing the number of flavors and fixing the number of 
quarks brings in some attractive channels. 

3-baryon potential (3 flavors}^ 

In this case the SU(3) decomposition is 

8 ® 8 ® 8 = 64 © (35 © 35) 2 © 27 6 © (10 © T0) 4 © 8 8 © 1 2 (24) 

and one can show that it is sufficient to work with the flavor content uuudddsss 
since all representations on the right hand side of the above formula occur in this 
flavor sector. 

The 3-baryon case is the most demanding technically. There are many Dirac 
index configurations to be considered and for some of them the dimension of the 
problem is quite large: with Dirac index structure 111223344 for example there are 
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Table 2. List of 3-baryon channels with anomalous dimensions greater than or equal to 3-/^. Only 
operators present at tree level are listed. 



Dirac structure 


48tt 2 7a 


SU(3) representation 


111223344 


48 


8 




44 


1,8 




42 


1,8 




36 


8,10,10 


111333224 


44 


1,8 


111222334 


48 


8 




44 


1 


111133442 


42 


1,8 


111122334 


44 


1 



originally 14130 operators and this number is reduced to 1518 after imposing the 
9-quark generalization of the identities (| 15[) . After performing the diagonalization 
in all Dirac sectors we found that in the 3-baryon case, although the vast majority 
of the eigenvalues are < 37jv, there are also some attractive channels. We enumerate 
these cases in Table [2] (only operators appearing at tree level in the OPE are listed). 

Our calculations allow us to determine the leading UV behavior of the total 3- 
baryon potential V%b- However, to draw any conclusion concerning the UV limit of 
the genuine 3-body potential V$b is difficult since comparing the results in Tables[T] 
and [2] we see that in some cases the 2-body forces are dominant, moreover there are 
also cases where there is cancellation between two different 2-body contributions. 

4. Results from lattice QCD and comparisons with OPE 
predictions 

In lattice QCD simulations with the spatial volume L 3 , the NBS wave function can 
be extracted from the 4-point correlation functionplEI which is given, for example, 
in the case of the TV TV system at t > to by 

Q a p(x,t-to) = (0\T{N a (r,t)N p {r + x,t)}J NN (t o )\0) (25) 

= Yl Mkn,s 1 , S2 )^ s ^(x)e^J t ^ + --- (26) 

n,si ,S2 

^A(k , Sl ,s 2 )i P k a ^(x)e^o (t ~ t o\ t^t , (27) 

with the matrix element A(k n , Sx, S2) = i n (NN, k n , Si, s 2 \J^nn(0)), where J^NN{to) 
is some source operator which can create states \NN, fe„, sx, S2)in, an d the ellipses 
represent inelastic contributions from intermediate states other than NN. In our 
numerical simulations, we have mainly used a wall source with the Coulomb gauge 
fixing only at t = io-^The large t — to necessary to achieve ground state saturation 
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in the last line causes some difficulties in numerical simulations. The signal-to- 
noise ratio for nucleon 4-pt functions behaves asymptotically for large t eJr^^Sl 

yjj^J ~ e _2 ( m « _3m "/ 2 ) t 5 where m^v and m, are the nucleon and pion masses, 

respectively. This problem becomes worse as we decrease the pion mass toward its 

physical value. Furthermore, as we increase the volume, the splitting between the 

ground state and the 1st excited state for the 2 nucleon system becomes smaller 

k 2 1 /2n\ 2 

as AE ~ = — . If L ~ 5/m and m N ~ 1 GeV, AE ~ 62MeV ~ 

niN rriN \ L J 

1/(3. 2fm), which requires t (AE)^ 1 ~ 3.2 fm for ground state saturation. The 
behavior of statistical noise in the above, however, makes the signals very poor at 
such large t for the two nucleon system. 

Recently, an improved extraction of the NBS wave function has been proposed to 
overcome the above difficulties! 2 ^ We first define the normalized 4-point correlation 
function as 

R(xA) = g(x,t)/(e- m " t ) 2 ~ ]T A(k n , Sl ,s 2 ) ip k -^(x)e- Aw kJ (28) 

n.s±.S2 

where AWk = 2\J m^ + k 2 — 2mjv- Using an identity AWk = k 2 /win — 

(AWk) 2 1 '(4mA/-) and neglecting inelastic contributions, we arrive at the time- 
dependent Schrodinger-like equation 

{-Ho - | + ^} R(x,t) = J * V U(*,vWv,t) (29) 

~ V(x)R(x,t) +■■■ , (30) 

which shows that the same U (x, y) in eq. (U]) can be obtained from R(x, t), where in 
the last line the derivative expansion is employed. An advantage of this extraction is 
that ground state saturation is not required for R(x, t) to satisfy eq. (|2T)|) or eq. (|3T))) . 
For this method to work, t has to be large enough such that elastic contributions 
dominate R(x, t). 

4.1. Nucleon- Nucleon potentials 

We first show the parity-even NN potentials at the leading order in the derivative 
expansion in Fig.[TJ where the central potential for the spin-singlet sector V r / =1 (r) = 
V Q I=1 (r) — 3V r / =1 (r), the central potential for the spin-triplet sector V r / =0 (r) = 
V I=0 (r) + Vj =Q (r) and the tensor potential for the spin-triplet sector V^. =0 (r) are 
plotted as a function of r by blue, green and red symbols, respectively.^ The results 
in the left figure have been obtained using PACS-CS gauge configurations in the 
2+1 flavor QCD at ~ 701 MeV and the lattice spacing a ~ 0.091 fm on 32 3 x 64 
latticepS while those in the right figure have been obtained in quenched QCD at 
m-r ~ 731 MeV and a ~ 0.137 fm on a 32 4 lattice. 

The central potentials for both spin-singlet and spin-triplet sectors have repul- 
sive cores, which confirm eq. (|17|) with positive R, while the tensor potential is less 
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Fig. 1. Central and tensor potentials in 2+1 flavor QCD at m n = 701 MeV by the improved 
method (left) and in quenched QCD at = 731 MeV by the conventional method(right). 

singular at short distance, in accordance with the fact that short distance contri- 
butions do not exist at the leading order of PT.ESl Note however that the central 
potentials are finite at the origin in contrast to the divergent behavior predicted 
by the OPE analysis, probably due to the artifact of non-zero lattice spacing. It is 
therefore important to show the appearance of the expected divergent behavior in 
future numerical simulations, by taking the continuum limit. We also observe that 
the repulsive cores and the tensor potential become significantly enhanced in the 
2+1 flavor QCD. 

4.2. Baryon interactions in the flavor SU(3) limit 

The figures in Fig.[2]show the leading order BB potential from full QCD simulations 
in the flavor SU(3) limit at the pseudo-scalar meson mass M-ps = 469(l)MeV and 
a = 0.121(2)fmj^2l where the central potentials (red) in the spin-singlet sector are 
given for the 27, 8 S and 1 irreducible representations of SU(3), from left to right 
in the 1st row, while the central and tensor potentials in the spin-triplet sector are 
given for 10, 10 and 8 Q in the 2nd row. 

The OPE analysis in the previous section predicts that, while the central po- 
tentials in the 27, 10 and 10 representations have a repulsive core, attractive in- 
teractions dominate at short distances for the central potentials in the 1, 8 S and 
8 a representations. Indeed, the central potential in the singlet representation, V^ x \ 
shows attraction at short distance in Fig. O which leads to the bound H-dibaryon 
in this channel.^ The OPE analysis correctly predicts the attractive interaction at 
short distance in the singlet representation. On the other hand, while the repulsion 
at short distance becomes weaker in the 8 a representation, it becomes strongest 
among all 6 representations in the 8 S representation. The OPE analysis does not 
reproduce these repulsions. 

This puzzle may be resolved by observing^ that there is no 6-quark operator in 
the 8 S channel in the nonrelativistic limit and therefore the corresponding nonper- 
turbative matrix elements must be small. Although the channel remains attractive 
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Fig. 2. Potentials of baryon-baryon S-wave interaction in the flavor 5(7(3) limit, labeled by the 
flavor irreducible representation. These are obtained at the pseudo-scalar meson mass of 469 McV 
by the improved method. 



at extremely short distances, in the range between 0.1 fm and 0.5 fm relevant for 
lattice simulations subleading operators, which are repulsive, dominate. Note that 
a valence quark model with gluon exchange^^l can reproduce the correct pattern 
of the short distance behavior of potentials observed in lattice QCD simulations 
including the attraction in the singlet channel. 



4.3. Three-nucleon force 

The three-nucleon force (3NF) has been calculated in the linear setup with p = 0, 
where three nucleons are aligned linearly with equal spacing ri = \t\/2, in 2-flavor 
QCD at a ~ 0.156 fm and m w ~ 1.13 GeV on a 16 3 x32 latticePln the left of Fig.[3j 
we give the 3NF V^tr, p) at p = as a function of which are obtained by the 
improved method in eqs. (|29ll30|) for the 3NF^"^ where R(x, t) for two nucleons is 
replaced with R(r, p, t) for three nucleons. The results from the different sink times 
are consistent with each other within statistical errors. We observe an indication of 
repulsive 3NF at short distance, which is consistent with the prediction by the OPE 
analysis that 3NF is universally repulsive. As in the case of the 2NF, the potential is 
finite at the origin due to the non-zero lattice spacing in lattice QCD simulations. 
It is therefore important to further investigate the short distance behavior as a 
function of the lattice spacing a. 
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Fig. 3. (Left) The 3NF Vzn{T> p) at p = as a function of T2 = [t*|/2, obtained by the improved 
method at (t — to)/a = 7.5(green), 8(red) and 8.5(blue) at m„- ~ 1.13 GeV and a ~ 0.156 fm. 
(Right) The spin-singlet central potential in the 2-flavor QCD with m w ~ 1.13 GeV as a function 
of r at a ~ 0.108 fm (red) and 0.156 fm (blue). 



5. Conclusions and discussions 

In this article, we have reviewed the recent activities to determine the short distance 
behavior of potentials among baryons, which are defined through the NBS wave 
function in the HAL QCD method. We can determine the short distance behaviors 
of the NBS wave functions and thus the corresponding potentials, employing the 
OPE (operator product expansion) and the RG analysis analytically in perturbation 
theory of QCD, thanks to its asymptotic freedom. 

The most notable predictions from the OPE analysis are (1) an existence of the 
repulsive cores for the two nucleon potentials with the use of a simple effective model 
for the ratio of the matrix elements, (2) an appearance of the attractive core, instead 
of the repulsive one, for the flavor singlet channel in the 3-flavor QCD, and (3) an 
existence of universal repulsion at short distance for the 3NF without relying on any 
non-perturbative inputs. These 3 results more or less agree with the behaviors of 
the corresponding potentials obtained in lattice QCD simulations, except the lattice 
potentials are always finite at the origin due to the lattice artifacts. There are, 
however, some disagreements between the OPE analysis and lattice QCD results, 
which are probably caused by the effects of the neglected subleading terms in the 
OPE analysis at not asymptotically small distances. 

The above success of the OPE analysis suggests that the repulsion among 
baryons at short distance may be explained by the combination of the Pauli suppres- 
sion principle among quarks and the one-gluon exchange between quarks^. Indeed 
in the OPE analysis, the tendency of repulsions become stronger for more numbers 
of quarks, while it becomes weaker for a larger number of valence flavors with the 
total number of quarks fixed. Since the OPE analysis also predicts forms of the 
repulsive/ attractive cores as a function of distances, it is very important and in- 

b The valence quark model with gluon exchange mentioned before relies on these properties .I21IHI 
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teresting in the future to confirm them in lattice QCD simulations by taking the 
continuum limit. The right of Fig. [3] shows our very preliminary result for a com- 
parison of the spin-singlet central potential between two lattice spacings, a ~ 0.108 
fm(red) and 0.156 fm(blue). We indeed observe an expected behavior that the re- 
pulsive core increases as we decrease the lattice spacing. More data and analysis in 
this direction will be required to determine the r dependence of the repulsive cores. 
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